1 Preparació de les dades (PRA1)


1.1 Presentació de les dades i establiment de l’objectiu analític


La cardiotocografia és una tècnica utilitzada per fer un monitoratge del ritme cardíac del fetus durant l’embaràs i el part. Es tracta d’una operació de monitoratge que, a part de la pulsació, pren informació de diverses variables per tal de detectar anomalies en el ritme del batec del cor, i determinar l’existència d’anomalies és per tant una tasca complicada de determinar de forma visual en temps real. Un dels aparells més utilitzats per fer aquest monitoratge és l’Omniview-SisPorto, un sistema automatitzat que realitza una anàlisi de la cardiotocografia i activa senyals d’alerta visual i sonora quan la combinació de mesures formen un patró considerat sospitós o patològic.

Malgrat el seu bon funcionament, s’han donat alguns casos en què patrons anòmals del ritme cardíac fetal han passat desapercebuts, i l’empresa Omniview ens ha contactat amb l’encàrrec de millorar els mètodes de classificació que fa el sistema SisPorto. L’objectiu és crear un model de classificació que millori els índexs actuals d’especificitat del sistema, de forma que es minimitzi el nombre de falsos negatius en la detecció d’anomalies del ritme cardíac fetal.

Per començar la nostra tasca, partim del dataset “Fetal cardiotocography data”, disponible a https://www.kaggle.com/datasets/akshat0007/fetalhr.

Es tracta d’un joc de dades amb 2129 observacions de 40 variables, que inclou dades capturades amb l’aparell SisPorto, però també altres variables recollides de forma manual per l’expert sanitari que realitzà la cardiotocografia en cada observació. Les dades inclouen també un diagnòstic que conclou, per part de l’expert mèdic, si hi ha un patró anòmal o potencialment patològic en cada observació. Aquest diagnòstic es reflecteix en dues de les variables del joc de dades. La primera, anomenada ‘Class’, és una escala de factors codificats de l’1 al 10, anant de normal a sospitós. La segona variable, també adient per aplicar algorismes de classificació, anomenada ‘NSP’, classifica les observacions en tres classes: 1=Normal, 2=Sospitós i 3=Patològic.


Presentem la totalitat de les variables que apareixen al joc de dades. En primer lloc tenim les variables que identifiquen les observacions.

-FileName (string): nom del fitxer de la cardiotocografia

-Data (string): data (mm/dd/aaaa) de l’exploració

-SegFile (string): identificador únic per a cada mesura.


Seguidament, tenim dues variables que determinen la duració (en segons) del monitoratge:

-b (numèrica): moment d’inici

-e (numèrica): moment de finalització


Les variables enregistrades pel sistema SisPorto:

-LB (numèrica): nivell (ritme cardíac) de base

-AC (numèrica): acceleracions

-FM (numèrica): moviments fetals

-UC (numèrica): contraccions uterines

-ASTV (numèrica): percentatge de temps amb variabilitat anormal de curta durada

-mSTV (numèrica): valor mitjà de variabilitat de curta durada

-ALTV (numèrica): percentatge de temps amb variabilitat anormal de llarga durada

-mLTV (numèrica): valor mitjà de variabilitat de llarga durada


Altres variables referents al ritme cardíac:

-LB (numèrica): nivell de base determinat per un expert

-DL (numèrica): desacceleracions lleugeres

-DS (numèrica): desacceleracions severes

-DP (numèrica): desacceleracions perllongades

-DR (numèrica): desacceleracions repetitives


Més enllà, el registre del ritme cardíac fetal genera un histograma, les dades del qual s’inclouen al dataset en les següents variables:

-Width (numèrica): amplada de l’histograma

-Min (numèrica): mínim de freqüència de l’histograma

-Max (numèrica): màxim de freqüència de l’histograma

-Nmax (numèrica): nombre de pics de l’histograma

-Nzeros (numèrica): nombre de zeros de l’histograma

-Mode (numèrica): mode

-Mean (numèrica): mitjana

-Median (numèrica): mediana

-Variance (numèrica): variància

-Tendency (categòrica): -1=asimètrica a l’esquerra, 0=simètrica, 1=asimètrica a la dreta


Finalment, hi ha un seguit de variables amb valoracions diagnòstiques de classificació:

-A (binària): son calmat

-B (binària): son REM

-C (binària): vetlla calmada

-D (binària): vetlla activa

-AD (binària): patró d’estrès

-DE (binària): patró de desacceleració

-LD (binària): patró sever de desacceleració

-FS (binària): patró sinusoïdal pla (estat patològic)

-SUSP (binària): patró sospitós

-CLASS (categòrica): codi de 0 a 1 per classes A a SUSP

-NSP (categòrica): 1=Normal, 2=Sospitós, 3=Patològic



1.2 Verificació i preparació de les dades


Carreguem les dades:

CTG <- read.csv("CTG.csv")


Mostrem les primeres línies del dataset:

as.data.frame(head(CTG))
##       FileName      Date     SegFile   b    e LBE  LB AC FM UC ASTV MSTV ALTV
## 1 Variab10.txt 12/1/1996 CTG0001.txt 240  357 120 120  0  0  0   73  0.5   43
## 2   Fmcs_1.txt  5/3/1996 CTG0002.txt   5  632 132 132  4  0  4   17  2.1    0
## 3   Fmcs_1.txt  5/3/1996 CTG0003.txt 177  779 133 133  2  0  5   16  2.1    0
## 4   Fmcs_1.txt  5/3/1996 CTG0004.txt 411 1192 134 134  2  0  6   16  2.4    0
## 5   Fmcs_1.txt  5/3/1996 CTG0005.txt 533 1147 132 132  4  0  5   16  2.4    0
## 6   Fmcs_2.txt  5/3/1996 CTG0006.txt   0  953 134 134  1  0 10   26  5.9    0
##   MLTV DL DS DP DR Width Min Max Nmax Nzeros Mode Mean Median Variance Tendency
## 1  2.4  0  0  0  0    64  62 126    2      0  120  137    121       73        1
## 2 10.4  2  0  0  0   130  68 198    6      1  141  136    140       12        0
## 3 13.4  2  0  0  0   130  68 198    5      1  141  135    138       13        0
## 4 23.0  2  0  0  0   117  53 170   11      0  137  134    137       13        1
## 5 19.9  0  0  0  0   117  53 170    9      0  137  136    138       11        1
## 6  0.0  9  0  2  0   150  50 200    5      3   76  107    107      170        0
##   A B C D E AD DE LD FS SUSP CLASS NSP
## 1 0 0 0 0 0  0  0  0  1    0     9   2
## 2 0 0 0 0 0  1  0  0  0    0     6   1
## 3 0 0 0 0 0  1  0  0  0    0     6   1
## 4 0 0 0 0 0  1  0  0  0    0     6   1
## 5 0 1 0 0 0  0  0  0  0    0     2   1
## 6 0 0 0 0 0  0  0  1  0    0     8   3


Verifiquem les dimensions del dataset:

dim(CTG)
## [1] 2129   40



Mostrem un resum de les variables:

summary(CTG)
##    FileName             Date             SegFile                b         
##  Length:2129        Length:2129        Length:2129        Min.   :   0.0  
##  Class :character   Class :character   Class :character   1st Qu.:  55.0  
##  Mode  :character   Mode  :character   Mode  :character   Median : 538.0  
##                                                           Mean   : 878.4  
##                                                           3rd Qu.:1521.0  
##                                                           Max.   :3296.0  
##                                                           NA's   :3       
##        e             LBE              LB              AC        
##  Min.   : 287   Min.   :106.0   Min.   :106.0   Min.   : 0.000  
##  1st Qu.:1009   1st Qu.:126.0   1st Qu.:126.0   1st Qu.: 0.000  
##  Median :1241   Median :133.0   Median :133.0   Median : 1.000  
##  Mean   :1703   Mean   :133.3   Mean   :133.3   Mean   : 2.722  
##  3rd Qu.:2435   3rd Qu.:140.0   3rd Qu.:140.0   3rd Qu.: 4.000  
##  Max.   :3599   Max.   :160.0   Max.   :160.0   Max.   :26.000  
##  NA's   :3      NA's   :3       NA's   :3       NA's   :3       
##        FM                UC              ASTV            MSTV      
##  Min.   :  0.000   Min.   : 0.000   Min.   :12.00   Min.   :0.200  
##  1st Qu.:  0.000   1st Qu.: 1.000   1st Qu.:32.00   1st Qu.:0.700  
##  Median :  0.000   Median : 3.000   Median :49.00   Median :1.200  
##  Mean   :  7.503   Mean   : 3.669   Mean   :47.01   Mean   :1.335  
##  3rd Qu.:  2.000   3rd Qu.: 5.000   3rd Qu.:61.00   3rd Qu.:1.700  
##  Max.   :564.000   Max.   :23.000   Max.   :87.00   Max.   :7.000  
##  NA's   :2         NA's   :2        NA's   :2       NA's   :2      
##       ALTV             MLTV              DL               DS          
##  Min.   : 0.000   Min.   : 0.000   Min.   : 0.000   Min.   :0.000000  
##  1st Qu.: 0.000   1st Qu.: 4.600   1st Qu.: 0.000   1st Qu.:0.000000  
##  Median : 0.000   Median : 7.400   Median : 0.000   Median :0.000000  
##  Mean   : 9.885   Mean   : 8.208   Mean   : 1.576   Mean   :0.003759  
##  3rd Qu.:11.000   3rd Qu.:10.800   3rd Qu.: 3.000   3rd Qu.:0.000000  
##  Max.   :91.000   Max.   :50.700   Max.   :16.000   Max.   :1.000000  
##  NA's   :2        NA's   :2        NA's   :1        NA's   :1         
##        DP               DR        Width             Min              Max     
##  Min.   :0.0000   Min.   :0   Min.   :  3.00   Min.   : 50.00   Min.   :122  
##  1st Qu.:0.0000   1st Qu.:0   1st Qu.: 37.00   1st Qu.: 67.00   1st Qu.:152  
##  Median :0.0000   Median :0   Median : 67.50   Median : 93.00   Median :162  
##  Mean   :0.1278   Mean   :0   Mean   : 70.45   Mean   : 93.58   Mean   :164  
##  3rd Qu.:0.0000   3rd Qu.:0   3rd Qu.:100.00   3rd Qu.:120.00   3rd Qu.:174  
##  Max.   :4.0000   Max.   :0   Max.   :180.00   Max.   :159.00   Max.   :238  
##  NA's   :1        NA's   :1   NA's   :3        NA's   :3        NA's   :3    
##       Nmax            Nzeros             Mode            Mean      
##  Min.   : 0.000   Min.   : 0.0000   Min.   : 60.0   Min.   : 73.0  
##  1st Qu.: 2.000   1st Qu.: 0.0000   1st Qu.:129.0   1st Qu.:125.0  
##  Median : 3.000   Median : 0.0000   Median :139.0   Median :136.0  
##  Mean   : 4.068   Mean   : 0.3236   Mean   :137.5   Mean   :134.6  
##  3rd Qu.: 6.000   3rd Qu.: 0.0000   3rd Qu.:148.0   3rd Qu.:145.0  
##  Max.   :18.000   Max.   :10.0000   Max.   :187.0   Max.   :182.0  
##  NA's   :3        NA's   :3         NA's   :3       NA's   :3      
##      Median         Variance         Tendency             A         
##  Min.   : 77.0   Min.   :  0.00   Min.   :-1.0000   Min.   :0.0000  
##  1st Qu.:129.0   1st Qu.:  2.00   1st Qu.: 0.0000   1st Qu.:0.0000  
##  Median :139.0   Median :  7.00   Median : 0.0000   Median :0.0000  
##  Mean   :138.1   Mean   : 18.81   Mean   : 0.3203   Mean   :0.1806  
##  3rd Qu.:148.0   3rd Qu.: 24.00   3rd Qu.: 1.0000   3rd Qu.:0.0000  
##  Max.   :186.0   Max.   :269.00   Max.   : 1.0000   Max.   :1.0000  
##  NA's   :3       NA's   :3        NA's   :3         NA's   :3       
##        B                C                 D                E          
##  Min.   :0.0000   Min.   :0.00000   Min.   :0.0000   Min.   :0.00000  
##  1st Qu.:0.0000   1st Qu.:0.00000   1st Qu.:0.0000   1st Qu.:0.00000  
##  Median :0.0000   Median :0.00000   Median :0.0000   Median :0.00000  
##  Mean   :0.2723   Mean   :0.02493   Mean   :0.0381   Mean   :0.03387  
##  3rd Qu.:1.0000   3rd Qu.:0.00000   3rd Qu.:0.0000   3rd Qu.:0.00000  
##  Max.   :1.0000   Max.   :1.00000   Max.   :1.0000   Max.   :1.00000  
##  NA's   :3        NA's   :3         NA's   :3        NA's   :3        
##        AD               DE               LD                FS         
##  Min.   :0.0000   Min.   :0.0000   Min.   :0.00000   Min.   :0.00000  
##  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.00000   1st Qu.:0.00000  
##  Median :0.0000   Median :0.0000   Median :0.00000   Median :0.00000  
##  Mean   :0.1562   Mean   :0.1185   Mean   :0.05033   Mean   :0.03246  
##  3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:0.00000   3rd Qu.:0.00000  
##  Max.   :1.0000   Max.   :1.0000   Max.   :1.00000   Max.   :1.00000  
##  NA's   :3        NA's   :3        NA's   :3         NA's   :3        
##       SUSP             CLASS            NSP       
##  Min.   :0.00000   Min.   : 1.00   Min.   :1.000  
##  1st Qu.:0.00000   1st Qu.: 2.00   1st Qu.:1.000  
##  Median :0.00000   Median : 4.00   Median :1.000  
##  Mean   :0.09266   Mean   : 4.51   Mean   :1.304  
##  3rd Qu.:0.00000   3rd Qu.: 7.00   3rd Qu.:1.000  
##  Max.   :1.00000   Max.   :10.00   Max.   :3.000  
##  NA's   :3         NA's   :3       NA's   :3




Si seleccionem les files amb valors NA veiem que es tracta de les últimes tres files del dataset, on la majoria de camps estan buits:

CTG[rowSums(is.na(CTG)) > 0,]
##      FileName Date SegFile  b  e LBE LB AC  FM UC ASTV MSTV ALTV MLTV DL DS DP
## 2127                       NA NA  NA NA NA  NA NA   NA   NA   NA   NA NA NA NA
## 2128                       NA NA  NA NA NA  NA NA   NA   NA   NA   NA  0  0  0
## 2129                       NA NA  NA NA NA 564 23   87    7   91 50.7 16  1  4
##      DR Width Min Max Nmax Nzeros Mode Mean Median Variance Tendency  A  B  C
## 2127 NA    NA  NA  NA   NA     NA   NA   NA     NA       NA       NA NA NA NA
## 2128  0    NA  NA  NA   NA     NA   NA   NA     NA       NA       NA NA NA NA
## 2129  0    NA  NA  NA   NA     NA   NA   NA     NA       NA       NA NA NA NA
##       D  E AD DE LD FS SUSP CLASS NSP
## 2127 NA NA NA NA NA NA   NA    NA  NA
## 2128 NA NA NA NA NA NA   NA    NA  NA
## 2129 NA NA NA NA NA NA   NA    NA  NA



Així doncs, eliminem aquests darrers tres registres sense perdre informació valuosa:

CTG2 <- na.omit(CTG)
dim(CTG2)
## [1] 2126   40



Comprovem que amb aquesta operació, el dataset queda net de valors buits:

sum(is.na(CTG2))
## [1] 0


Seguidament, observem que hi ha dos variables que no ens semblen d’especial rellevància de cara a la tasca. Es tracta de les variables ‘FileName’ i ‘Date’. Considerem, d’una banda, que els registres queden correctament identificats si mantenim la variable ‘SegFile’, i d’altra banda, la data del registre no hauria de tenir cap influència en la classificació (assumint que la precisió de l’instrument de mesura no es veu afectat per la data de la mesura). Així doncs, eliminem aquestes dues variables del dataset:

CTG2 <- subset(CTG2, select = -c(FileName, Date))


Seguidament, obtenint informació suplementària sobre el joc de dades al repositori de UCI, ens adonem que les variables ‘AC’, ‘FM’, ‘UC’, ‘DL’, ‘DS’ i ‘DP’ apareixen al dataset com a nombres absoluts d’ocurrències. Per tal que aquestes variables es puguin comparar entre registres, cal que estiguin expressades en nombre d’ocurrències per segon. Com que cada registre té una durada diferent, cal fer aquesta operació registre a registre. Per fer-ho, utilitzem les variables ‘b’ i ‘e’ (moment d’inici i de finalització) per determinar la durada del registre, i creem noves variables que expressen ocurrències per segon:

CTG2$ACps <- round(CTG2$AC / (CTG2$e - CTG2$b), digits = 3)
CTG2$FMps <- round(CTG2$FM / (CTG2$e - CTG2$b), digits = 3)
CTG2$UCps <- round(CTG2$UC / (CTG2$e - CTG2$b), digits = 3)
CTG2$DLps <- round(CTG2$DL / (CTG2$e - CTG2$b), digits = 3)
CTG2$DSps <- round(CTG2$DS / (CTG2$e - CTG2$b), digits = 3)
CTG2$DPps <- round(CTG2$DP / (CTG2$e - CTG2$b), digits = 3)


Mostrem un nou resum per seguir explorant les dades:

str(CTG2)
## 'data.frame':    2126 obs. of  44 variables:
##  $ SegFile : chr  "CTG0001.txt" "CTG0002.txt" "CTG0003.txt" "CTG0004.txt" ...
##  $ b       : int  240 5 177 411 533 0 240 62 120 181 ...
##  $ e       : int  357 632 779 1192 1147 953 953 679 779 1192 ...
##  $ LBE     : int  120 132 133 134 132 134 134 122 122 122 ...
##  $ LB      : int  120 132 133 134 132 134 134 122 122 122 ...
##  $ AC      : int  0 4 2 2 4 1 1 0 0 0 ...
##  $ FM      : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ UC      : int  0 4 5 6 5 10 9 0 1 3 ...
##  $ ASTV    : int  73 17 16 16 16 26 29 83 84 86 ...
##  $ MSTV    : num  0.5 2.1 2.1 2.4 2.4 5.9 6.3 0.5 0.5 0.3 ...
##  $ ALTV    : int  43 0 0 0 0 0 0 6 5 6 ...
##  $ MLTV    : num  2.4 10.4 13.4 23 19.9 0 0 15.6 13.6 10.6 ...
##  $ DL      : int  0 2 2 2 0 9 6 0 0 0 ...
##  $ DS      : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ DP      : int  0 0 0 0 0 2 2 0 0 0 ...
##  $ DR      : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Width   : int  64 130 130 117 117 150 150 68 68 68 ...
##  $ Min     : int  62 68 68 53 53 50 50 62 62 62 ...
##  $ Max     : int  126 198 198 170 170 200 200 130 130 130 ...
##  $ Nmax    : int  2 6 5 11 9 5 6 0 0 1 ...
##  $ Nzeros  : int  0 1 1 0 0 3 3 0 0 0 ...
##  $ Mode    : int  120 141 141 137 137 76 71 122 122 122 ...
##  $ Mean    : int  137 136 135 134 136 107 107 122 122 122 ...
##  $ Median  : int  121 140 138 137 138 107 106 123 123 123 ...
##  $ Variance: int  73 12 13 13 11 170 215 3 3 1 ...
##  $ Tendency: int  1 0 0 1 1 0 0 1 1 1 ...
##  $ A       : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ B       : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ C       : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ D       : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ E       : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ AD      : int  0 1 1 1 0 0 0 0 0 0 ...
##  $ DE      : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ LD      : int  0 0 0 0 0 1 1 0 0 0 ...
##  $ FS      : int  1 0 0 0 0 0 0 1 1 1 ...
##  $ SUSP    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ CLASS   : int  9 6 6 6 2 8 8 9 9 9 ...
##  $ NSP     : int  2 1 1 1 1 3 3 3 3 3 ...
##  $ ACps    : num  0 0.006 0.003 0.003 0.007 0.001 0.001 0 0 0 ...
##  $ FMps    : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ UCps    : num  0 0.006 0.008 0.008 0.008 0.01 0.013 0 0.002 0.003 ...
##  $ DLps    : num  0 0.003 0.003 0.003 0 0.009 0.008 0 0 0 ...
##  $ DSps    : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ DPps    : num  0 0 0 0 0 0.002 0.003 0 0 0 ...


Tenim la impressió que pot ser que no hi hagi diferència entre la línia de base mesurada pel sistema SisPorto i la línia de base establerta per l’expert:

all(CTG2$LBE == CTG2$LB)
## [1] TRUE


Efectivament, les dues variables són idèntiques, i per tant en podem eliminar una de les dues:

CTG2 <- subset(CTG2, select = -LBE)


Ens adonem també que tots els valors de la variable ‘DR’ són 0, i per tant també eliminem aquesta variable:

summary(CTG2$DR)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       0       0       0       0
CTG3 <- subset(CTG2, select = -DR)


1.2.1 Detecció d’outliers


Hem analitzat una per una les variables numèriques que considerem interessants de cara a una possible tasca de classificació. A continuació mostrem algunes de les variables on hem detectat valors aïllats:


La variable ‘AC’ (acceleracions), mostra un gran nombre de valors aïllats:

boxplot.stats(CTG3$AC)$out
##  [1] 12 11 13 12 11 15 11 13 12 17 13 11 17 11 18 12 11 11 13 16 14 14 12 11 13
## [26] 12 14 12 14 12 12 11 12 12 19 19 13 17 13 13 15 12 14 26 18 12 17 13 14 13
## [51] 17 14 14 14 11 11 12 13 11 11 14 13 14 11 15 13 15 14 13 12 15 14 16 13 12
## [76] 17 17 16 16 11 12 11 21
boxplot(CTG3$AC, main = "AC")


Si mostrem el valor més allunyat, veiem que la resta de variables no mostren valors estranys, i per tant la totalitat del registre sembla tenir coherència. Es podria tractar d’un error només en el valor de la variable ‘AC’. Fixant-nos-hi més en detall, veiem que la duració del registre és de més de 3000 segons, una duració molt llarga, que fa que un nombre total gran d’ocurrències de la variable ‘AC’ no sigui tan estrany:

CTG3[CTG3$AC > 25,]
##          SegFile   b    e  LB AC FM UC ASTV MSTV ALTV MLTV DL DS DP Width Min
## 1648 CTG1650.txt 357 3591 130 26  3 13   51  1.6    0  4.2 12  0  0   109  77
##      Max Nmax Nzeros Mode Mean Median Variance Tendency A B C D E AD DE LD FS
## 1648 186    4      0  144  138    143       38        0 0 0 0 0 0  1  0  0  0
##      SUSP CLASS NSP  ACps  FMps  UCps  DLps DSps DPps
## 1648    0     6   1 0.008 0.001 0.004 0.004    0    0


Per comparar-ho, fem la mateixa operació de detecció d’outliers, però aquest cop utilitzem la variable ‘ACps’, que mostra ocurrències per segon en lloc de nombre absolut d’ocurrències:

boxplot.stats(CTG3$ACps)$out
##  [1] 0.017 0.016 0.019 0.016 0.016 0.016 0.017 0.016 0.018 0.017 0.018 0.016
## [13] 0.017 0.016


boxplot(CTG3$ACps, main = "ACps")

Veiem que els valors extrems són menys, i també menys allunyats. Considerem que aquests valors allunyats no són errors de registre, i els mantenim en el dataset sense cap correcció. Realitzem la resta de l’anàlisi de detecció d’outliers utilitzant les variables que registren ocurrències per segon:


CTG3[CTG3$ACps > 0.016,]
##          SegFile    b    e  LB AC FM UC ASTV MSTV ALTV MLTV DL DS DP Width Min
## 182  CTG0182.txt  194  948 138 13  0  4   35  5.3    0  4.5  0  0  0   148  52
## 530  CTG0530.txt  829 1192 142  7 31  0   32  2.3    0  0.0  0  0  0   144  56
## 631  CTG0631.txt  136  941 134 14  2  3   48  2.2    0  0.0  0  0  0   120  50
## 1095 CTG1095.txt 1548 2114 122 10  0  1   22  2.5    0  2.2  0  0  0    52 100
## 1097 CTG1097.txt 1352 1871 123  9  0  1   24  2.2    0  1.7  0  0  0    52 100
## 1249 CTG1250.txt 2337 2897 112 10  0  4   25  1.4    0  0.6  0  0  0    35 109
## 1860 CTG1862.txt  828 1648 138 14  0  3   51  0.9    0  0.2  0  0  0    49 122
##      Max Nmax Nzeros Mode Mean Median Variance Tendency A B C D E AD DE LD FS
## 182  200   11      2  146  157    161       72        1 0 0 0 1 0  0  0  0  0
## 530  200   10      0  170  158    162       37        1 0 0 0 1 0  0  0  0  0
## 631  170    5      0  160  150    155       28        1 0 0 0 1 0  0  0  0  0
## 1095 152    1      0  136  132    134        6        0 0 1 0 0 0  0  0  0  0
## 1097 152    1      0  136  133    135        5        0 0 1 0 0 0  0  0  0  0
## 1249 144    0      0  116  120    119        6       -1 0 1 0 0 0  0  0  0  0
## 1860 171    4      0  147  149    150        5        0 0 1 0 0 0  0  0  0  0
##      SUSP CLASS NSP  ACps  FMps  UCps DLps DSps DPps
## 182     0     4   1 0.017 0.000 0.005    0    0    0
## 530     0     4   1 0.019 0.085 0.000    0    0    0
## 631     0     4   1 0.017 0.002 0.004    0    0    0
## 1095    0     2   1 0.018 0.000 0.002    0    0    0
## 1097    0     2   1 0.017 0.000 0.002    0    0    0
## 1249    0     2   1 0.018 0.000 0.007    0    0    0
## 1860    0     2   1 0.017 0.000 0.004    0    0    0


Analitzem la variable ‘FMps’ (moviments fetals per segon):

boxplot.stats(CTG3$FMps)$out
##   [1] 0.072 0.222 0.408 0.380 0.441 0.383 0.451 0.469 0.340 0.425 0.334 0.135
##  [13] 0.099 0.108 0.112 0.089 0.103 0.085 0.109 0.079 0.065 0.055 0.058 0.047
##  [25] 0.038 0.012 0.018 0.020 0.009 0.010 0.008 0.008 0.028 0.026 0.107 0.009
##  [37] 0.008 0.013 0.016 0.029 0.050 0.053 0.058 0.009 0.011 0.008 0.013 0.013
##  [49] 0.016 0.008 0.015 0.008 0.010 0.022 0.026 0.015 0.015 0.010 0.010 0.021
##  [61] 0.008 0.016 0.018 0.013 0.017 0.016 0.019 0.008 0.009 0.009 0.010 0.009
##  [73] 0.008 0.009 0.010 0.009 0.013 0.013 0.018 0.015 0.012 0.019 0.025 0.012
##  [85] 0.020 0.018 0.020 0.021 0.020 0.014 0.021 0.022 0.025 0.022 0.024 0.028
##  [97] 0.016 0.015 0.019 0.024 0.009 0.019 0.011 0.016 0.019 0.021 0.010 0.019
## [109] 0.008 0.019 0.022 0.009 0.009 0.009 0.009 0.009 0.011 0.013 0.013 0.016
## [121] 0.019 0.013 0.008 0.008 0.009 0.012 0.009 0.017 0.023 0.023 0.020 0.020
## [133] 0.024 0.035 0.054 0.030 0.048 0.088 0.030 0.054 0.043 0.009 0.054 0.011
## [145] 0.058 0.052 0.085 0.091 0.026 0.033 0.092 0.084 0.115 0.084 0.015 0.017
## [157] 0.013 0.014 0.014 0.015 0.035 0.033 0.029 0.041 0.024 0.040 0.027 0.029
## [169] 0.031 0.026 0.063 0.060 0.050 0.085 0.071 0.060 0.029 0.008 0.008 0.008
## [181] 0.020 0.010 0.027 0.017 0.018 0.010 0.008 0.010 0.009 0.014 0.010 0.013
## [193] 0.025 0.009 0.013 0.017 0.017 0.021 0.025 0.018 0.013 0.017 0.010 0.008
## [205] 0.306 0.298 0.139 0.189 0.014 0.157 0.235 0.408 0.360 0.455 0.443 0.470
## [217] 0.477 0.446 0.481 0.369 0.335 0.430 0.346 0.323 0.375 0.346 0.353 0.010
## [229] 0.021 0.017 0.045 0.016 0.011 0.008 0.017 0.017 0.020 0.030 0.012 0.043
## [241] 0.035 0.040 0.048 0.012 0.014 0.013 0.011 0.008 0.010 0.010 0.012 0.012
## [253] 0.012 0.014 0.008 0.010 0.012 0.009 0.011 0.010 0.014 0.011 0.010 0.010
## [265] 0.010 0.010 0.010 0.011 0.010 0.016 0.010 0.008 0.009 0.008 0.011 0.011
## [277] 0.017 0.019 0.019 0.025 0.028 0.048 0.043 0.032 0.040 0.051 0.036 0.037
## [289] 0.050 0.035 0.049 0.049 0.041 0.050 0.048 0.054 0.052 0.008 0.008 0.015
## [301] 0.017 0.008 0.009 0.009 0.009 0.010 0.009


boxplot(CTG3$FMps, main = "FMps")


hist(CTG3$FMps, breaks = sqrt(nrow(CTG3)))

En el cas de la variable FMps, hi ha en un principi una quantitat extrema d’outliers. L’anàlisi visual mostra que en realitat no es tracta d’outliers, sinó que estem davant d’una variable amb una distribució molt asimètrica. La variable ‘FM’ mesura ‘moviments fetals’, i sembla lògic que hi hagi una gran variabilitat, així com una gran quantitat de valors baixos o zero, ja que quan el fetus dorm o està en estat de vetlla calmada, els moviments són pocs. En tot cas, creem una nova variable amb una transformació logarítmica d’aquesta variable, en cas que ens pugui resultar útil en models posteriors:

library(ggplot2)
CTG3$FMps_log <- log(CTG3$FMps)

ggplot(mapping = aes(CTG3$FMps_log)) + geom_density()
## Warning: Removed 1311 rows containing non-finite values (stat_density).


Analitzem la variable ‘UCps’(contraccions uterines per segon):

boxplot.stats(CTG3$UCps)$out
## [1] 0.015


boxplot(CTG3$UCps, main = "UCps")


Detectem un únic valor allunyat, però no tenim evidència per considerar-lo com a erroni, i per tant el mantenim en el dataset:

CTG3[CTG3$UCps > 0.014,]
##          SegFile   b    e  LB AC FM UC ASTV MSTV ALTV MLTV DL DS DP Width Min
## 1165 CTG1165.txt 910 1379 131  5  0  7   26  1.5    0    3  0  0  0    61 109
##      Max Nmax Nzeros Mode Mean Median Variance Tendency A B C D E AD DE LD FS
## 1165 170    2      1  155  151    154       11        1 0 0 0 1 0  0  0  0  0
##      SUSP CLASS NSP  ACps FMps  UCps DLps DSps DPps FMps_log
## 1165    0     4   1 0.011    0 0.015    0    0    0     -Inf



Analitzem la variable ‘ASTV’ (percentatge de temps amb variabilitat anormal de curta durada), i no hi trobem outliers:

boxplot.stats(CTG3$ASTV)$out
## integer(0)



Analitzem la variable ‘ALTV’ (percentatge de temps amb variabilitat anormal de llarga durada), i ens trobem amb la mateixa situació que amb la variable ‘FM’, una distribució molt asimètrica, amb molts valors allunyats, però que no es poden considerar com a registres erronis. No hi fem cap modificació:

boxplot.stats(CTG3$ALTV)$out
##   [1] 43 79 72 71 40 69 54 53 38 29 37 29 67 68 75 74 30 49 39 31 32 34 38 31 58
##  [26] 39 46 57 75 29 33 40 51 37 34 34 54 58 52 32 38 59 78 84 71 59 52 62 39 45
##  [51] 30 45 35 67 72 45 58 32 41 61 56 84 56 84 44 39 71 71 61 67 81 42 57 56 61
##  [76] 77 67 88 91 78 84 47 41 59 55 61 49 49 82 86 29 67 78 84 60 68 30 32 41 28
## [101] 38 29 49 35 62 68 70 47 48 51 31 43 50 47 71 48 39 39 38 31 31 34 28 29 38
## [126] 32 37 32 37 36 28 28 31 40 72 85 64 58 66 91 90 54 70 77 58 81 29 40 38 54
## [151] 73 43 50 42 53 64 41 34 33 34 37 44 62 29 69 72 91 90 84 42 91 66 40 34 32
## [176] 30 32 38 32 45 32 32 32 31 33 39 39 44 55 41 33 31 32 32 35 53 44 31 33 57
## [201] 50 64 59 67 62 71 74 59 48 45 62 59 50 61 39 43 56 62 53 44 48 61 54 37 58
## [226] 40 30 37 42 28 47 48 28 43 38 28 39 32 61 48 35 34 45 49 42 33 32 32 31 33
## [251] 35 32 31 28 33 36 49 52 46 41 37 54 60 62 59 62 56 49 38 32 45 30 44 30 46
## [276] 37 40 54 65 46 34 44 46 55 58 59 64 60 61 73 71 67 63 52 58 35 28 73 70 59
## [301] 44 34 29 40 36 33 48 36 36
boxplot(CTG3$ALTV, main = "ALTV")



Pel que fa a la variable ‘DS’, ens trobem amb una situació diferent. Els possibles valors són 0 i 1, i dels 2126 casos, només 7 es classifiquen com a ‘1’. La variable mesura ‘desacceleracions severes’, i és per tant normal que només una petita minoria dels casos presentin aquest patró:

CTG3[CTG3$DS != 0,]
##          SegFile    b    e  LB AC FM UC ASTV MSTV ALTV MLTV DL DS DP Width Min
## 1489 CTG1491.txt 1532 2655 132  2  0  9   31  1.4    0 11.5  0  1  1   102  61
## 1490 CTG1492.txt 1728 2655 132  0  0  6   32  1.3    0 13.6  0  1  1    91  60
## 1792 CTG1794.txt 1839 2894 121  0  1  4   66  2.1    0  6.4 11  1  0   105  55
## 1793 CTG1795.txt 1922 2894 121  0  1  3   67  2.1    0  0.0 11  1  0   102  55
## 1794 CTG1796.txt 1922 2771 121  0  1  4   66  2.1    0  0.0 10  1  0   102  55
## 1795 CTG1797.txt 2018 2892 121  0  1  3   68  2.1    0  0.0  9  1  0   102  55
## 1796 CTG1798.txt 2153 2892 121  0  0  3   70  1.9    0  0.0  7  1  0   102  55
##      Max Nmax Nzeros Mode Mean Median Variance Tendency A B C D E AD DE LD FS
## 1489 163    5      0   99  121    129       94        1 0 0 0 0 0  1  0  0  0
## 1490 151    1      1   99  116    125       72        1 0 0 0 0 0  0  0  1  0
## 1792 160    7      0   67   85     92      109       -1 0 0 0 0 0  0  0  1  0
## 1793 157    4      1   67   81     87       89       -1 0 0 0 0 0  0  0  1  0
## 1794 157    5      1   67   83     90       98       -1 0 0 0 0 0  0  0  1  0
## 1795 157    3      1   67   79     82       83       -1 0 0 0 0 0  0  0  1  0
## 1796 157    6      2   67   76     79       68       -1 0 0 0 0 0  0  0  1  0
##      SUSP CLASS NSP  ACps  FMps  UCps  DLps  DSps  DPps  FMps_log
## 1489    0     6   1 0.002 0.000 0.008 0.000 0.001 0.001      -Inf
## 1490    0     8   3 0.000 0.000 0.006 0.000 0.001 0.001      -Inf
## 1792    0     8   3 0.000 0.001 0.004 0.010 0.001 0.000 -6.907755
## 1793    0     8   3 0.000 0.001 0.003 0.011 0.001 0.000 -6.907755
## 1794    0     8   3 0.000 0.001 0.005 0.012 0.001 0.000 -6.907755
## 1795    0     8   3 0.000 0.001 0.003 0.010 0.001 0.000 -6.907755
## 1796    0     8   3 0.000 0.000 0.004 0.009 0.001 0.000      -Inf


1.2.2 Discretització


Procedim a discretitzar diverses de les variables. En primer lloc, afegirem etiquetes a la variable ‘DS’ que acabem d’analitzar, creant una nova variable amb els codis ‘absent’ i ‘present’:

library(arules)
## Loading required package: Matrix
## 
## Attaching package: 'arules'
## The following objects are masked from 'package:base':
## 
##     abbreviate, write
CTG3$DS_d <- discretize(CTG3$DS, method = "fixed", breaks = c(-Inf, 1, Inf),labels = c("absent", "present"))
table(CTG3$DS_d)
## 
##  absent present 
##    2119       7


Seguidament discretitzarem totes les variables que mesuren ocurrències per segon, creant en cada cas tres nivells: ‘baix’, ‘mitjà’ i ‘alt’. Farem la discretització utilitzant el mètode k-means:

CTG3$AC_d <- discretize(CTG3$ACps, method = "cluster", labels = c("low", "medium", "high"))
CTG3$FM_d <- discretize(CTG3$FMps, method = "cluster", labels = c("low", "medium", "high"))
CTG3$UC_d <- discretize(CTG3$UCps, method = "cluster", labels = c("low", "medium", "high"))
CTG3$DL_d <- discretize(CTG3$DLps, method = "cluster", labels = c("low", "medium", "high"))
CTG3$DP_d <- discretize(CTG3$DPps, method = "cluster", labels = c("low", "medium", "high"))
summary(CTG3$AC_d)
##    low medium   high 
##   1037    660    429
summary(CTG3$FM_d)
##    low medium   high 
##   2037     61     28
summary(CTG3$UC_d)
##    low medium   high 
##    450    616   1060
summary(CTG3$DL_d)
##    low medium   high 
##   1627    404     95
summary(CTG3$DP_d)
##    low medium   high 
##   1948     70    108



També crearem una nova variable ‘NSP2’ a partir de la variable ‘NSP’, on ajuntarem les categories ‘sospitós’ i ‘patològic’. Amb això obtenim una variable target classificatòria binària, i podem més endavant triar si els nostres models de classificació tenen 2 o 3 classes a identificar:

summary(as.factor(CTG3$NSP))
##    1    2    3 
## 1655  295  176
CTG3$NSP2 <- as.factor(CTG3$NSP)
levels(CTG3$NSP2) <- c('Normal', 'Suspect', 'Pathologic')
levels(CTG3$NSP2)[levels(CTG3$NSP2)=="Pathologic"] <-"Suspect"
summary(CTG3$NSP2)
##  Normal Suspect 
##    1655     471


Finalment, i abans de donar per acabat el treball de neteja de les dades, ens assegurem que totes les variables categòriques o binàries estiguin factoritzades:

factcols <- c('A', 'B', 'C', 'D', 'AD', 'DE', 'LD', 'FS', 'SUSP', 'CLASS', 'NSP')
CTG3[factcols] <- lapply(CTG3[factcols], factor)
sapply(CTG3, class)
##     SegFile           b           e          LB          AC          FM 
## "character"   "integer"   "integer"   "integer"   "integer"   "integer" 
##          UC        ASTV        MSTV        ALTV        MLTV          DL 
##   "integer"   "integer"   "numeric"   "integer"   "numeric"   "integer" 
##          DS          DP       Width         Min         Max        Nmax 
##   "integer"   "integer"   "integer"   "integer"   "integer"   "integer" 
##      Nzeros        Mode        Mean      Median    Variance    Tendency 
##   "integer"   "integer"   "integer"   "integer"   "integer"   "integer" 
##           A           B           C           D           E          AD 
##    "factor"    "factor"    "factor"    "factor"   "integer"    "factor" 
##          DE          LD          FS        SUSP       CLASS         NSP 
##    "factor"    "factor"    "factor"    "factor"    "factor"    "factor" 
##        ACps        FMps        UCps        DLps        DSps        DPps 
##   "numeric"   "numeric"   "numeric"   "numeric"   "numeric"   "numeric" 
##    FMps_log        DS_d        AC_d        FM_d        UC_d        DL_d 
##   "numeric"    "factor"    "factor"    "factor"    "factor"    "factor" 
##        DP_d        NSP2 
##    "factor"    "factor"


Creem un nou dataframe amb totes les transformacions realitzades:

ctg1 <- CTG3




write.csv(ctg1, "ctg1.csv", row.names = FALSE)


ctgp2 <- read.csv("ctg1.csv")



2 Respostes PRA 2


2.1 Exercici 1


Com vam comentar a la pràctica anterior, l’objectiu del nostre estudi és millorar la capacitat de classificació de models construïts amb les dades proporcionades per sistema SisPorto. El nostre joc de dades conté una gran quantitat de dades, tant numèriques com categòriques, però seguint el nostre objectiu, seleccionem les variables numèriques que genera els sistema SisPorto per aplicar-hi un model no supervisat. Prescindim de DSps (desacceleracions severes), ja que entenem que quan es detecta una desacceleració severa, això per si sol ja indica un cas patològic. Ens interessa en canvi veure com les altres variables, que capturen anomalies més lleus o subtils, es combinen dintre de cada grup diagnòstic.

num_df <- ctgp2[c('ACps', 'FMps', 'UCps', 'DLps', 'DPps', 'ASTV', 'ALTV')]


Normalitzem les dades:

#definim la funció de normalització
mm_norm <- function(x) {
    (x - min(x)) / (max(x) - min(x))
}

num_norm <- as.data.frame(lapply(num_df, mm_norm))



if (!require('cluster')) install.packages('cluster')
## Loading required package: cluster
library(cluster)

Per tal de començar a trobar el nombre ideal de clústers en un model no supervisat, iterem la funció kmeans per agrupacions d’entre 2 i 10 clústers:

d <- daisy(num_norm) 
results <- rep(0, 10)
for (i in c(2,3,4,5,6,7,8,9,10))
{
  fit           <- kmeans(num_norm, i)
  y_cluster     <- fit$cluster
  sk            <- silhouette(y_cluster, d)
  results[i] <- mean(sk[,3])
}



Prenent com a mesura l’índex de silhouette, el nombre òptim de clústers seria 3, que precisament coincideix amb els tres grups que tenim definits segons la variable NSP (‘normal’, ‘sospitós’ i ‘patològic’):

plot(2:10,results[2:10],type="o",col="blue",pch=0,xlab="Number of clusters",ylab="Silhouette")




Si, com a alternativa, fem el càlcul de la suma de quadrats de distàncies respecte a centroides, també per a entre 2 i 10 agrupacions, el gràfic resultant és més difícil d’interpretar que l’anterior, però sembla que hi pugui haver un ‘colze’ també a l’alçada de 3 clústers:

results <- rep(0, 10)
for (i in c(2,3,4,5,6,7,8,9,10))
{
  fit           <- kmeans(num_norm, i)
  results[i] <- fit$tot.withinss
}
plot(2:10,results[2:10],type="o",col="blue",pch=0,xlab="Nombre de clústers",ylab="tot.tot.withinss")


Pel que fa als índexs ‘Calinski-Harabasz’ (‘ch’) i silueta mitjana (‘asw’), aquests ens dónen respectivament un nombre òptim de clústers de 2 i 3:

if (!require('fpc')) install.packages('fpc'); library('fpc')
## Loading required package: fpc
fit_ch  <- kmeansruns(num_norm, krange = 1:10, criterion = "ch") 
fit_asw <- kmeansruns(num_norm, krange = 1:10, criterion = "asw") 
fit_ch$bestk
## [1] 2
fit_asw$bestk
## [1] 3


A partir d’aquest moment, per avaluar l’algorisme de classificació, tenim dues opcions. Podem comparar la classificació amb 3 clústers comparant-la amb la variable NSP, que determina 3 nivells, o també podem crear un model amb 2 clústers, i comparar-lo amb la variable NSP2 que havíem creat anteriorment, i que ajuntava els nivells ‘sospitós’ i ‘patològic’ en un de sol.

Veiem una primera comparació entre les variables ALTV i ASTV amb 3 clústers:

ctg3clusters <- kmeans(num_norm, 3)

plot(num_norm[c(6,7)], col=ctg3clusters$cluster, main="Classificació k-means")

Com podem veure, si comparem aquest gràfic amb la classificació real, els tres clústers creats per l’algorisme no corresponen de forma gaire fidel als 3 grups diagnòstics. Hi ha una mena d’inversió dels grups. El grup definit pels nivells baixos d’ambdues variables en la classificació real (en negre al gràfic de sota) es veu dividit en 2 subgrups en la classificació amb k-means, i els dos altres grups reals, en la zona de valors progressivament més alts, es veuen compactats en un sol grup.

plot(num_norm[c(6,7)], col=as.factor(ctgp2$NSP), main="Classificació real")



Provem ara de comparar la classificació de k-means de 2 clústers amb la variable NSP2, que utilitza només 2 etiquetes diagnòstiques.

ctg2clusters <- kmeans(num_norm, 2)

plot(num_norm[c(6,7)], col=ctg2clusters$cluster, main="Classificació k-means")

plot(num_norm[c(6,7)], col=as.factor(ctgp2$NSP2), main="Classificació real")

A la comparació dels gràfics veiem que, a grans trets, la classificació sembla coincidir, si més no pels valors extrems de les dues variables ALTV i ASTV. Però si ens hi fixem, en la zona intermèdia, la pertinença a un grup o l’altre difereix notablement entre les dues gràfiques. Mentre que a la classificació real l’encavalcament és més gradual, i fins i tot hi ha casos pertanyents al primer grup que s’enfilen fins els valors més alts de les dues variables, a la classificació amb k-means els dos grups es veuen dividits sense pràcticament cap encavalcament.

Realitzem el mateix tipus de comparació per les variables FMps i ACps:

plot(num_norm[c(1,2)], col=ctg3clusters$cluster, main="Classificació k-means")

plot(num_norm[c(1,2)], col=as.factor(ctgp2$NSP), main="Classificació real")
legend("topright",                   
       legend = c("normal", "sospitós", "patològic"),
       col = 1:3, pch = 1)



plot(num_norm[c(1,2)], col=ctg2clusters$cluster, main="Classificació k-means")

plot(num_norm[c(1,2)], col=as.factor(ctgp2$NSP2), main="Classificació real")
legend("topright",                   
       legend = c("normal", "sospitós"),
       col = 1:2, pch = 1)

Tant en 2 com en 3 clústers, la divisió que proposa k-means difereix visiblement de les etiquetes diagnòstiques.

Fem una altra comparació, aquest cop només per 3 clústers, amb les variables ALTV i ACps:

plot(num_norm[c(1,7)], col=ctg3clusters$cluster, main="Classificació k-means")

plot(num_norm[c(1,7)], col=as.factor(ctgp2$NSP), main="Classificació real")
legend("topright",                   
       legend = c("normal", "sospitós", "patològic"),
       col = 1:3, pch = 1)

A la classificació real veiem com la categoria ‘normal’ es distribueix tot al llarg dels valors mitjans baixos de ACps, i desapareix per complet en els valors alts d’ALTV. En canvi, l’algorisme k-means situa dos grups diferents al llarg de la variable ACps, i també situa un tercer grup ocupant tots els valors de la variable ALTV que tenen 0 com a valors d’ACps.


Fem una última comparació entre les variables ALTV i FMps:

plot(num_norm[c(2,7)], col=ctg3clusters$cluster, main="Classificació k-means")

plot(num_norm[c(2,7)], col=as.factor(ctgp2$NSP), main="Classificació real")
legend("topright",                   
       legend = c("normal", "sospitós", "patològic"),
       col = 1:3, pch = 1)

Veiem un efecte semblant a la comparació anterior. A la classificació real, els valors baixos, propers a 0, de la variable FMps, mostren registres distribuïts a través de la variable ALTV de forma que podem identificar els tres grups diferents. El grup ‘normal’ té, o bé valor 0 de FMps, o bé valor 0 de ALTV. En canvi, a la classificació amb k-means, trobem un grup amb valors 0 de FMps, i els altres dos grups amb valors 0 d’ALTV. Es tracta doncs d’una classificació força diferent en el cas de la relació entre aquestes dues variables.




2.2 Exercici 2


Ara generarem un segon model de classificació no supervisada, però aquest cop utilitzarem la fórmula de correlació de Pearson absoluta com a mètrica per mesurar la distància.


library(amap)
ctg3clusters_abpearson <- Kmeans(num_norm, 3, method = "abspearson")
## Warning: did not converge in 10 iterations


Comparem el nou model, amb 3 clústers, amb la classificació real, per les variables ALTV i ASTV:

plot(num_norm[c(6,7)], col=ctg3clusters_abpearson$cluster, main="Classificació k-means")

plot(num_norm[c(6,7)], col=as.factor(ctgp2$NSP), main="Classificació real")

Com en el model anterior, els 3 clústers que genera el model no s’acaben d’ajustar al tres grups diagnòstics.

Ara compararem, per les variables ALTV i FMps, el primer model que ja teníem, el model amb absolute Pearson, i la classificació real:

plot(num_norm[c(2,7)], col=ctg3clusters$cluster, main="Classificació k-means")

plot(num_norm[c(2,7)], col=ctg3clusters_abpearson$cluster, main="Classificació k-means/abspearson")

plot(num_norm[c(2,7)], col=as.factor(ctgp2$NSP), main="Classificació real")
legend("topright",                   
       legend = c("normal", "sospitós", "patològic"),
       col = 1:3, pch = 1)

Als gràfics veiem com, malgrat que els dos models k-means no són idèntics, tendeixen a dividir els clústers de forma similar, i cap dels dos s’aproxima a la classificació real.

Fem el mateix tipus de comparació per les variables ASTV i DPps, on veurem el mateix efecte:

plot(num_norm[c(5,6)], col=ctg3clusters_abpearson$cluster, main="Classificació k-means/abspearson")

plot(num_norm[c(5,6)], col=ctg3clusters$cluster, main="Classificació k-means")

plot(num_norm[c(5,6)], col=as.factor(ctgp2$NSP), main="Classificació real")
legend("topright",                   
       legend = c("normal", "sospitós", "patològic"),
       col = 1:3, pch = 1)



ctg2clusters_abpearson <- Kmeans(num_norm, 2, method = "abspearson")
## Warning: did not converge in 10 iterations



Ara fem una última comparació per les mateixes variables ASTV i DPps, però aquest cop comparem per 2 clústers:

plot(num_norm[c(5,6)], col=ctg2clusters$cluster, main="Classificació k-means")

plot(num_norm[c(5,6)], col=ctg2clusters_abpearson$cluster, main="Classificació k-means/abspearson")

plot(num_norm[c(5,6)], col=as.factor(ctgp2$NSP2), main="Classificació real")
legend("topright",                   
       legend = c("normal", "sospitós"),
       col = 1:2, pch = 1)

El primer algorisme k-means talla horitzontalment el gràfic: divideix clarament els grups segons els valors de la variable ASTV. Amb la mètrica abspearson, l’algorisme classifica en un clúster els valors superiors de ASTV, i també aquells que tenen valors alts d’ambdues variables. En aquest sentit, la divisió de clústers s’assembla més a la classificació real, tot i que no classifica bé els valors intermedis.






2.3 Exercici 3


if (!require('dbscan')) install.packages('dbscan'); library('dbscan')
## Loading required package: dbscan
## 
## Attaching package: 'dbscan'
## The following object is masked from 'package:fpc':
## 
##     dbscan


Fem una primera prova de l’algorisme DBSCAN amb un valor arbitrari d’eps = 1:

tryeps <- dbscan(num_norm, eps = 1)
tryeps
## DBSCAN clustering for 2126 objects.
## Parameters: eps = 1, minPts = 5
## The clustering contains 1 cluster(s) and 0 noise points.
## 
##    1 
## 2126 
## 
## Available fields: cluster, eps, minPts

L’algorisme genera un sol clúster. Caldrà cercar un valor eps diferent.


Generem un gràfic de distància k que ens pugui donar pistes del valor òptim d’eps:

kNNdistplot(num_norm, k = 5)
abline(h=.25, col = "red", lty=2)

Amb la línia vermella hem indicat la zona on es podria trobar el valor d’eps que ens interessa, que seria entre 0.2 i 0.3.

Després de diferents intents entorn a aquests valors, trobem que amb el valor eps=0.29 obtenim 3 clústers i només 7 registres cauen fora de la classificació:

ctg_dbs <- dbscan(num_norm, eps = .29)

ctg_dbs
## DBSCAN clustering for 2126 objects.
## Parameters: eps = 0.29, minPts = 5
## The clustering contains 3 cluster(s) and 7 noise points.
## 
##    0    1    2    3 
##    7 2088   23    8 
## 
## Available fields: cluster, eps, minPts



Ara utilitzarem l’algorisme OPTICS, i guardarem la variable que resulta d’extreure el dbscan per tal d’avaluar el model:

opti <- optics(num_norm)
opt <- extractDBSCAN(opti, eps_cl = .29)
opt
## OPTICS ordering/clustering for 2126 objects.
## Parameters: minPts = 5, eps = 0.742812272170372, eps_cl = 0.29, xi = NA
## The clustering contains 3 cluster(s) and 7 noise points.
## 
##    0    1    2    3 
##    7 2088    8   23 
## 
## Available fields: order, reachdist, coredist, predecessor, minPts, eps,
##                   eps_cl, xi, cluster



Visualitzem els resultats de la classificació. Recordem que, també en la nostra classificació real, tenim un grup molt nombrós (diagnòstic ‘normal’) i dos grups molt menys nombrosos (diagnòstics ‘sospitós’ i ‘patològic’):

summary(as.factor(ctg1$NSP))
##    1    2    3 
## 1655  295  176

Veiem que en els 3 clústers generats per l’algorisme optics aquesta diferència encara és més gran (només 8 i 23 registres en el 2n i 3r grup respectivament). Per aquest motiu, els gràfics no seran fàcils d’interpretar.

plot(opt)



Veiem una comparativa entre la classificació amb dbscan i la classificació real, per a totes les variables 2 a 2:

pairs(num_norm, col= opt$cluster)



pairs(num_norm, col= ctgp2$NSP)


Es fa difícil observar les diferències. Vegem doncs una comparació només per les variables ALTV i FMps:

plot(num_norm[c(2,7)], col=opt$cluster, main="Classificació DBSCAN")

plot(num_norm[c(2,7)], col=as.factor(ctgp2$NSP), main="Classificació real")
legend("topright",                   
       legend = c("normal", "sospitós", "patològic"),
       col = 1:3, pch = 1)


hullplot(num_norm[c(2,7)], opt)

Veiem en els gràfics que, comparant la classificació dbscan amb la classificació real, els resultats no semblen més gaire acurats respecte als models k-means que hem generat anteriorment.


Obtenim, amb l’índex de silhouette, una mètrica d’avaluació del model generat amb dbscan:

noise <- opt$cluster==0
clusters <- opt$cluster[!noise]
d <- dist(num_norm[!noise, 1:7])
silh <- silhouette(clusters, d)
plot(silh, border=NA, col=sort(clusters), main="")

Si tenim en compte que els valors de l’índex de silhouette se situen entre -1 i 1, el valor que obtenim de 0.34 indica una qualitat d’agrupament acceptable.

En definitiva, els tres models comparats fins ara (kmeans, kmeans amb mètrica de Pearson, i dbscan) generen grups que s’emmirallen d’alguna manera amb la classificació real segons l’etiqueta diagnòstica NSP, sobretot per que fa als valors extrems dels parells de variables, i també en el sentit que agrupen 3 clústers de mides similars (1 grup nombrós i 2 grups minoritaris). Malgrat això, cap dels tres models encerta a separar de forma gaire acurada els registres amb combinacions de valors intermedis en les diferents variables.




2.4 Exercici 4


Per la creació d’un model d’arbre de decisió, tenint en compte que l’aplicació i interpretació del model seran més fàcils i eficients si treballem amb variables categòriques, tenim la possibilitat de fer una nova selecció de variables. Tanmateix, per tal de seguir responent a l’objectiu inicial de recerca, ens interessa seguir incloent les variables que el sistema SisPorto mesura (recordem que no totes les variables de les quals disposem estan mesurades amb SisPorto). Així doncs, proposem fer algunes proves amb dos grups de variables. El primer inclourà la versió discretitzada de les variables mesurades amb SisPorto, a més de les variables diagnòstiques ‘AD’ (patró d’estrès) i ‘DE’ (patró de desacceleració), i el segon afegirà a aquesta mateixa selecció les variables diagnòstiques ‘A’ (son calmat), ‘B’ (son REM), ‘C’ (vetlla calmada), i ‘D’ (vetlla activa):

cat_ctg <- ctg1[c('AD', 'DE', 'AC_d', 'FM_d', 'UC_d', 'DL_d', 'DP_d', 'NSP2')]

cat_ctg_ext <- ctg1[c('A', 'B', 'C', 'D','AD', 'DE', 'AC_d', 'FM_d', 'UC_d', 'DL_d', 'DP_d', 'NSP2')]


Creem conjunts train i test per a les dues mostres:

set.seed(666)
y <- cat_ctg[,8] 
X <- cat_ctg[,1:7]
set.seed(666)
y2 <- cat_ctg_ext[,12] 
X2 <- cat_ctg_ext[,1:11]



split_prop <- 3 
indexes = sample(1:nrow(cat_ctg), size=floor(((split_prop-1)/split_prop)*nrow(cat_ctg)))
trainX<-X[indexes,]
trainy<-y[indexes]
testX<-X[-indexes,]
testy<-y[-indexes]
split_prop <- 3 
indexes = sample(1:nrow(cat_ctg_ext), size=floor(((split_prop-1)/split_prop)*nrow(cat_ctg_ext)))
trainX2<-X2[indexes,]
trainy2<-y2[indexes]
testX2<-X2[-indexes,]
testy2<-y2[-indexes]



Comprovem que la divisió de conjunts sigui equilibrada:

sum(summary(trainy))
## [1] 1417
sum(summary(trainX$AD))
## [1] 1417
sum(summary(testy))
## [1] 709
sum(summary(testX$DE))
## [1] 709
sum(summary(trainy2))
## [1] 1417
sum(summary(trainX2$A))
## [1] 1417
sum(summary(testy2))
## [1] 709
sum(summary(testX2$B))
## [1] 709



Mostrem les regles del primer arbre de decisió:

model1 <- C50::C5.0(trainX, trainy, rules=TRUE)
summary(model1)
## 
## Call:
## C5.0.default(x = trainX, y = trainy, rules = TRUE)
## 
## 
## C5.0 [Release 2.07 GPL Edition]      Tue Jun 14 23:21:50 2022
## -------------------------------
## 
## Class specified by attribute `outcome'
## 
## Read 1417 cases (8 attributes) from undefined.data
## 
## Rules:
## 
## Rule 1: (713/18, lift 1.2)
##  AC_d in {medium, high}
##  DP_d in {low, medium}
##  ->  class Normal  [0.973]
## 
## Rule 2: (136/6, lift 1.2)
##  DE = 1
##  AC_d = low
##  UC_d in {medium, high}
##  DP_d in {low, medium}
##  ->  class Normal  [0.949]
## 
## Rule 3: (1015/98, lift 1.1)
##  UC_d in {medium, high}
##  DP_d = low
##  ->  class Normal  [0.903]
## 
## Rule 4: (8, lift 4.3)
##  DE = 0
##  AC_d = low
##  DP_d = medium
##  ->  class Suspect  [0.900]
## 
## Rule 5: (68/8, lift 4.1)
##  DP_d = high
##  ->  class Suspect  [0.871]
## 
## Rule 6: (186/59, lift 3.2)
##  AC_d = low
##  UC_d = low
##  ->  class Suspect  [0.681]
## 
## Default class: Normal
## 
## 
## Evaluation on training data (1417 cases):
## 
##          Rules     
##    ----------------
##      No      Errors
## 
##       6  175(12.4%)   <<
## 
## 
##     (a)   (b)    <-classified as
##    ----  ----
##    1050    67    (a): class Normal
##     108   192    (b): class Suspect
## 
## 
##  Attribute usage:
## 
##   87.09% DP_d
##   85.53% UC_d
##   73.61% AC_d
##   10.16% DE
## 
## 
## Time: 0.0 secs

Comparem amb l’arbre de decisió per la segona mostra:

model2 <- C50::C5.0(trainX2, trainy2, rules=TRUE)
summary(model2)
## 
## Call:
## C5.0.default(x = trainX2, y = trainy2, rules = TRUE)
## 
## 
## C5.0 [Release 2.07 GPL Edition]      Tue Jun 14 23:21:50 2022
## -------------------------------
## 
## Class specified by attribute `outcome'
## 
## Read 1417 cases (12 attributes) from undefined.data
## 
## Rules:
## 
## Rule 1: (366, lift 1.3)
##  B = 1
##  ->  class Normal  [0.997]
## 
## Rule 2: (275/2, lift 1.3)
##  A = 1
##  ->  class Normal  [0.989]
## 
## Rule 3: (222/2, lift 1.3)
##  AD = 1
##  ->  class Normal  [0.987]
## 
## Rule 4: (54, lift 1.3)
##  D = 1
##  ->  class Normal  [0.982]
## 
## Rule 5: (40, lift 1.3)
##  C = 1
##  ->  class Normal  [0.976]
## 
## Rule 6: (158/14, lift 1.2)
##  DE = 1
##  FM_d in {low, medium}
##  ->  class Normal  [0.906]
## 
## Rule 7: (300/1, lift 4.4)
##  A = 0
##  B = 0
##  C = 0
##  D = 0
##  AD = 0
##  DE = 0
##  ->  class Suspect  [0.993]
## 
## Rule 8: (2, lift 3.3)
##  DE = 1
##  FM_d = high
##  ->  class Suspect  [0.750]
## 
## Default class: Normal
## 
## 
## Evaluation on training data (1417 cases):
## 
##          Rules     
##    ----------------
##      No      Errors
## 
##       8   19( 1.3%)   <<
## 
## 
##     (a)   (b)    <-classified as
##    ----  ----
##    1097     1    (a): class Normal
##      18   301    (b): class Suspect
## 
## 
##  Attribute usage:
## 
##   47.00% B
##   40.58% A
##   36.84% AD
##   32.46% DE
##   24.98% D
##   23.99% C
##   11.29% FM_d
## 
## 
## Time: 0.0 secs

El primer arbre de decisió genera 5 regles, i obté una estimació d’error entorn un 11%. (* els percentatges d’errors i índexs de precisió que s’indiquen d’ara en endavant són aproximats, ja que varien lleugerament en cada execució de l’algorisme, i per tant també varien quan exportem el document de Rmd a Html).

Les variables més utilitzades són DP_d (desacceleracions perllongades), FM_d (moviments fetals), i AD (patró d’estrès), que és una variable diagnòstica.

Pel que fa a la segona mostra, amb 12 variables, l’arbre de decisió genera 7 regles, amb una estimació d’error de només entre un 1% i un 2%.

Mostrem gràficament els arbres dels dos models:



model1_plot <- C50::C5.0(trainX, trainy)
plot(model1_plot)

predicted_model1 <- predict( model1, testX, type="class" )
print(sprintf("La precisió de l'arbre de decisió es: %.4f %%",100*sum(predicted_model1 == testy) / length(predicted_model1)))
## [1] "La precisió de l'arbre de decisió es: 84.4852 %"
model2_plot <- C50::C5.0(trainX2, trainy2)
plot(model2_plot)



predicted_model2 <- predict( model2, testX2, type="class" )
print(sprintf("La precisió de l'arbre de decisió del segon model és: %.4f %%",100*sum(predicted_model2 == testy2) / length(predicted_model2)))
## [1] "La precisió de l'arbre de decisió del segon model és: 98.5896 %"



D’aquests resultats, ens crida l’atenció la precisió del segon model. Sembla que és un model gairebé perfecte, amb una estimació d’error d’un 1.5% i una precisió de 98%. Això ens fa necessàriament sospitar de l’adequació del model. De fet, el primer model, amb un 86% de precisió, tampoc es queda curt.

La diferència entre els dos models té a veure amb la inclusió de variables diagnòstiques sumades a les variables mesurades per SisPorto. En el primer model en tenim dues, i en el segon model en tenim fins a 6. Una variable diagnòstica és, al capdavall, una variable classificatòria, i per tant té una funció similar a la variable NSP que estem utilitzant com a referència per supervisar el model. Pot ser doncs que això sigui un error de plantejament: si prenem com a exemple la variable AD (patró d’estrès), el sistema que classifica cada registre com a 0 o 1 en AD és el mateix sistema que classifica cada registre com a 1, 2 o 3 en la variable NSP, i és per tant poc sorprenent que la variable AD sigui predictora de NSP. Pot ser que el model sigui tan precís perquè reflecteix un efecte de circularitat: equivaldria a dir que un pacient està deprimit perquè té un trastorn depressiu, que és el mateix que dir que un pacient té un trastorn depressiu perquè està deprimit, quan el que ens interessa és establir que un pacient té un trastorn depressiu perquè té una desregulació del son, apatia, pensaments negatius, etc.

Així doncs, ens sembla que seria més correcte generar un model on prescindim de totes les variables classificatòries, i utilitzem només les variables mesurades amb SisPorto. Creiem que només així podem respondre al nostre objectiu inicial d’estudi. A més, això ens permetrà comparar els models no supervisats amb els models supervisats, ja que estarem treballant bàsicament amb les mateixes variables.



Hi ha dues variables mesurades per SisPorto, que hem utilitzat en els nostres models no supervisats, que encara no havíem discretitzat. Són ASTV (percentatge de temps amb variabilitat anormal de curta durada) i ALTV (percentatge de temps amb variabilitat anormal de llarga durada). Fem ara aquesta discretització:

library(arules)
ctg_plus <- ctg1
ctg_plus$ASTV_d <- discretize(ctg_plus$ASTV, method = "cluster", labels = c("low", "medium", "high"))
ctg_plus$ALTV_d <- discretize(ctg_plus$ALTV, method = "cluster", labels = c("low", "medium", "high"))

summary(ctg_plus$ASTV_d)
##    low medium   high 
##    655    619    852
summary(ctg_plus$ALTV_d)
##    low medium   high 
##   1704    275    147


Creem ara una tercera mostra que inclou aquestes dues variables i descarta totes les variables diagnòstiques, excepte NSP2, que serà la nostra variable classificatòria:

cat_ctg3 <- ctg_plus[c('ASTV_d', 'ALTV_d','AC_d', 'FM_d', 'UC_d', 'DL_d', 'DP_d', 'NSP2')]



Fem els passos previs per generar l’arbre de decisió:

set.seed(666)
y3 <- cat_ctg3[,8] 
X3 <- cat_ctg3[,1:7]
split_prop <- 3 
indexes = sample(1:nrow(cat_ctg3), size=floor(((split_prop-1)/split_prop)*nrow(cat_ctg3)))
trainX3<-X3[indexes,]
trainy3<-y3[indexes]
testX3<-X3[-indexes,]
testy3<-y3[-indexes]
sum(summary(trainy3))
## [1] 1417
sum(summary(trainX3$ASTV_d))
## [1] 1417
sum(summary(testy3))
## [1] 709
sum(summary(testX3$ALTV_d))
## [1] 709



Mostrem les regles d’aquest tercer arbre de decisió:

model3 <- C50::C5.0(trainX3, trainy3, rules=TRUE)
summary(model3)
## 
## Call:
## C5.0.default(x = trainX3, y = trainy3, rules = TRUE)
## 
## 
## C5.0 [Release 2.07 GPL Edition]      Tue Jun 14 23:21:51 2022
## -------------------------------
## 
## Class specified by attribute `outcome'
## 
## Read 1417 cases (8 attributes) from undefined.data
## 
## Rules:
## 
## Rule 1: (854/49, lift 1.2)
##  ASTV_d in {low, medium}
##  ->  class Normal  [0.942]
## 
## Rule 2: (1066/77, lift 1.2)
##  ALTV_d = low
##  DP_d in {low, medium}
##  ->  class Normal  [0.927]
## 
## Rule 3: (68/8, lift 4.1)
##  DP_d = high
##  ->  class Suspect  [0.871]
## 
## Rule 4: (283/120, lift 2.7)
##  ALTV_d in {medium, high}
##  ->  class Suspect  [0.575]
## 
## Default class: Normal
## 
## 
## Evaluation on training data (1417 cases):
## 
##          Rules     
##    ----------------
##      No      Errors
## 
##       4  151(10.7%)   <<
## 
## 
##     (a)   (b)    <-classified as
##    ----  ----
##    1066    51    (a): class Normal
##     100   200    (b): class Suspect
## 
## 
##  Attribute usage:
## 
##   95.20% ALTV_d
##   80.03% DP_d
##   60.27% ASTV_d
## 
## 
## Time: 0.0 secs

L’arbre genera 4 regles de decisió, obté una estimació d’error d’un 10%, i utilitza tres variables: ALTV, DP i ASTV.

Mostrem gràficament l’arbre de decisió:

model3_plot <- C50::C5.0(trainX3, trainy3)
plot(model3_plot)



Veiem que aquest model obté una precisió del 87%:

predicted_model3 <- predict( model3, testX3, type="class" )
print(sprintf("La precisió de l'arbre de decisió es: %.4f %%",100*sum(predicted_model3 == testy3) / length(predicted_model3)))
## [1] "La precisió de l'arbre de decisió es: 87.7292 %"



Ara mostrem en detall una matriu de confusió on podem comparar la predicció del model amb les dades del conjunt test:

if(!require(gmodels)){
    install.packages('gmodels', repos='http://cran.us.r-project.org')
    library(gmodels)
}
## Loading required package: gmodels



CrossTable(testy3, predicted_model3,prop.chisq  = FALSE, prop.c = FALSE, prop.r =FALSE,dnn = c('Reality', 'Prediction'))
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  709 
## 
##  
##              | Prediction 
##      Reality |    Normal |   Suspect | Row Total | 
## -------------|-----------|-----------|-----------|
##       Normal |       517 |        21 |       538 | 
##              |     0.729 |     0.030 |           | 
## -------------|-----------|-----------|-----------|
##      Suspect |        66 |       105 |       171 | 
##              |     0.093 |     0.148 |           | 
## -------------|-----------|-----------|-----------|
## Column Total |       583 |       126 |       709 | 
## -------------|-----------|-----------|-----------|
## 
## 



Ara creem un segon model que inclogui adaptive boosting:

model3_ab <- C50::C5.0(trainX3, trainy3, trials = 100)



Veiem que amb aquesta modificació obtenim una petita millora (d’entorn un 1%) en la precisió del model:

predicted_model3_ab <- predict( model3_ab, testX3, type="class" )
print(sprintf("La precisió del model amb adaptive boosting: %.4f %%",100*sum(predicted_model3_ab == testy) / length(predicted_model3_ab)))
## [1] "La precisió del model amb adaptive boosting: 88.7165 %"



Mostrem la matriu de confusió pel model amb adaptive boosting:

CrossTable(testy3, predicted_model3_ab, prop.chisq  = FALSE, prop.c = FALSE, prop.r =FALSE,dnn = c('Reality', 'Prediction'))
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  709 
## 
##  
##              | Prediction 
##      Reality |    Normal |   Suspect | Row Total | 
## -------------|-----------|-----------|-----------|
##       Normal |       511 |        27 |       538 | 
##              |     0.721 |     0.038 |           | 
## -------------|-----------|-----------|-----------|
##      Suspect |        53 |       118 |       171 | 
##              |     0.075 |     0.166 |           | 
## -------------|-----------|-----------|-----------|
## Column Total |       564 |       145 |       709 | 
## -------------|-----------|-----------|-----------|
## 
## 

Si comparem les dues matrius de confusió, veiem que la millora de precisió es reflecteix prinicipalment en la detecció de veritables negatius (categoria ‘Suspect’). Respecta al primer model, el model amb adaptive boosting prediu més casos sospitosos que realment ho són. Al mateix temps, el nombre de casos classificats erròniament com a normals es redueix amb adaptive boosting. Així doncs, malgrat que l’augment de precisió és només d’entorn els 2 punts percentuals, el segon model respon de forma més efectiva als nostres objectius.





2.5 Exercici 5


Ara generem un model amb l’algorisme Random Forest que classifica els registres segons les mateixes variables que hem utilitzat fins ara. Dividim els conjunts d’entrenament i test:

set.seed(123)

samp <- sample(nrow(cat_ctg3), 0.6669 * nrow(cat_ctg3))

train <- cat_ctg3[samp, ]

test <- cat_ctg3[-samp, ]



I apliquem el model:

library(randomForest)
## randomForest 4.7-1
## Type rfNews() to see new features/changes/bug fixes.
## 
## Attaching package: 'randomForest'
## The following object is masked from 'package:ggplot2':
## 
##     margin
model_rf <- randomForest(NSP2~., data = train)
model_rf
## 
## Call:
##  randomForest(formula = NSP2 ~ ., data = train) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 2
## 
##         OOB estimate of  error rate: 9.74%
## Confusion matrix:
##         Normal Suspect class.error
## Normal    1092      23   0.0206278
## Suspect    115     187   0.3807947

Obtenim una estimació d’error del 11%, havent utilitzat 2 variables a cada bifurcació del model.


Amb aquest model obtenim una precisió del 88%, un nivell similar als models d’arbre de decisió anteriors (lleugerament per sota del model amb adaptive boosting):

prediction_rf <- predict(model_rf, newdata = test)
sum(prediction_rf==test$NSP2) / nrow(test)
## [1] 0.8984485


A la representació gràfica del model veiem com les variables que es consideren com a més importants en la predicció del model són ASTV, ALTV, AC i DP (el nombre ’MeanDecreaseGini indica quin percentatge decauria la precisió de la predicció si alguna d’aquestes variables no es tingués en compte):

varImpPlot(model_rf)



Tot seguit, mostrem les matrius de confusió dels models d’arbre decisió sense i amb boosting, i del model Random Forest:

library(caret)
## Loading required package: lattice



  confusionMatrix(data= predicted_model3, reference = testy3)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction Normal Suspect
##    Normal     517      66
##    Suspect     21     105
##                                           
##                Accuracy : 0.8773          
##                  95% CI : (0.8509, 0.9005)
##     No Information Rate : 0.7588          
##     P-Value [Acc > NIR] : 1.772e-15       
##                                           
##                   Kappa : 0.6317          
##                                           
##  Mcnemar's Test P-Value : 2.390e-06       
##                                           
##             Sensitivity : 0.9610          
##             Specificity : 0.6140          
##          Pos Pred Value : 0.8868          
##          Neg Pred Value : 0.8333          
##              Prevalence : 0.7588          
##          Detection Rate : 0.7292          
##    Detection Prevalence : 0.8223          
##       Balanced Accuracy : 0.7875          
##                                           
##        'Positive' Class : Normal          
## 




confusionMatrix(data=predicted_model3_ab, reference = testy3)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction Normal Suspect
##    Normal     511      53
##    Suspect     27     118
##                                           
##                Accuracy : 0.8872          
##                  95% CI : (0.8615, 0.9095)
##     No Information Rate : 0.7588          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.6749          
##                                           
##  Mcnemar's Test P-Value : 0.005189        
##                                           
##             Sensitivity : 0.9498          
##             Specificity : 0.6901          
##          Pos Pred Value : 0.9060          
##          Neg Pred Value : 0.8138          
##              Prevalence : 0.7588          
##          Detection Rate : 0.7207          
##    Detection Prevalence : 0.7955          
##       Balanced Accuracy : 0.8199          
##                                           
##        'Positive' Class : Normal          
## 



confusionMatrix(data=prediction_rf, reference = test$NSP2)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction Normal Suspect
##    Normal     534      66
##    Suspect      6     103
##                                           
##                Accuracy : 0.8984          
##                  95% CI : (0.8738, 0.9197)
##     No Information Rate : 0.7616          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.6815          
##                                           
##  Mcnemar's Test P-Value : 3.57e-12        
##                                           
##             Sensitivity : 0.9889          
##             Specificity : 0.6095          
##          Pos Pred Value : 0.8900          
##          Neg Pred Value : 0.9450          
##              Prevalence : 0.7616          
##          Detection Rate : 0.7532          
##    Detection Prevalence : 0.8463          
##       Balanced Accuracy : 0.7992          
##                                           
##        'Positive' Class : Normal          
## 

Ja havíem comparat els dos primers models, i vèiem que el model amb adaptive boosting obté major precisió, però sobretot obté una especificitat molt més gran (classifica correctament com a casos sospitosos els que realment ho són). El model Random Forest no millora els resultats del model amb adaptive boosting, tot i que la seva sensibilitat és molt alta.





2.6 Exercici 6


El nostre estudi tenia com a objectiu plantejar models amb una alta capacitat de classificació. Concretament, ens interessava optimitzar l’especificitat dels possibles models de selecció, ja que imaginàvem el cas en què en l’anàlisi d’algunes cardiotocografies, el sistema SisPorto havia diagnosticat com a normal alguns fetus amb patrons anòmals de ritme cardíac.

El nostre joc de dades tenia fins a 40 variables de diferents tipus, incloent les mesures del sistema SisPorto, però també altres variables amb càlculs derivats dels registres, i variables categòriques diagnòstiques.

Per tal de centrar-nos en l’objectiu d’anàlisi, i també per poder comparar els resultats dels diferents models, ens hem cenyit a les variables mesurades pel sistema SisPorto tant pels models no supervisats com pels models supervisats (hem discretitzat variables per adaptar-les als models supervisats).

Hem constatat que els models no supervisats no mostren, en el cas del nostre joc de dades, una gran capacitat de classificar correctament segons els 3 grups diagnòstics (variable NSP: ‘normal’, ‘sospitós’ i ‘patològic’). Val a dir però que els models no supervisats k-means sí que identifica un nombre òptim de clústers (2 o 3) que té relació amb aquesta mateixa variable classificatòria que nosaltres coneixem d’entrada. Quan, per avaluar els models no supervisats que hem generat (k-means, k-means amb mètrica absolute Pearson i DBSCAN), fem comparacions entre l’agrupació de l’algorisme i les dades reals per 2 variables, veiem repetidament que els models tenen dificultats per discriminar correctament els registres que mostren valors intermedis de les variables. Una altra dificultat, que es fa visible especialment en l’aplicació del model DBSCAN, és que els nostres 3 grups diagnòstics tenen mides molt desiguals: els casos patològics representen una part petita dels registres, i la majoria de casos són diagnosticats com a normals. El model DBSCAN exagera aquesta desigualtat, i classifica com a patològics una quantitat excessivament petita de registres. En definitiva, un model no supervisat, almenys entre els que hem estudiat, no seria una bona opció per detectar casos patològics de forma automàtica.

Pel que fa als models supervisats, hem trobat que, en general, aconseguíem uns índexs de precisió molt alts. La quantitat de registres i variables semblava formar una bona base per l’entrenament dels algorismes, tant pel que fa als arbres de decisió com l’algorisme Random Forest. En ambdós casos, no semblava tampoc que hi hagués un problema de sobreentrenament, ja que els algorismes seguien mostrant bons resultats sobre els conjunts test. Hem constatat que, entre els 3 models supervisats que hem generat, el model d’arbre de decisió amb adaptive boosting mostra el nivell d’especificitat més alt, i classifica erròniament casos patològics amb menor freqüència comparat amb els altres dos models. Tanmateix, si imaginem utilitzar aquest algorisme per detectar de forma automàtica els casos patològics, trobem que encara hi ha massa casos que es classifiquen com a normals que en realitat són patològics. I aquí rau el perill d’utilitzar aquests models en el nostre cas concret. Tot i que hem pogut generar models supervisats amb bons índexs de precisió, cal tenir present la rellevància de l’error que suposa diagnosticar com a normal una cardiotocografia d’un fetus que té un patró patològic de ritme cardíac. Això no significa necessàriament que els nostre models, o d’altres més òptims que es puguin generar, no siguin útils, sinó que cal seguir utilitzant-los com a suport o contrast a l’avaluació d’un expert que pugui tenir en compte altres variables més enllà de les mesurades per sistema automàtic SisPorto.